//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
///////////////////////////////////////////////////////////////////////
// template<class T> class SymMatrixSolver
// is a specialized class for a symmetric matrix eigen-value problems.
// Using Lapack Library
// usage (e.g., see OrthogonalTB::getDynamicMatrix)
//  SymMatrixSolver<double> dmsolver;
//  dmsolver.solve(mydatafirst, ndim);
///////////////////////////////////////////////////////////////////////
#ifndef SYMMETRICMATRIXSOLVER_H
#define SYMMETRICMATRIXSOLVER_H
#include "EigSolver/EigenSpaceDimension.h"
#include "EigSolver/LapackSymSolver.h"

template<class T>
class SymMatrixRef {
public:
  typedef T Type_t;
  SymMatrixRef(T* a):d_ptr(a) { }
  ~SymMatrixRef() { }
  T* begin() { return d_ptr;}

private:
  T* d_ptr;
};

template<class T>
class SymMatrixSolver {

public:

  typedef 
  LapackSymSolver<SymMatrixRef<T>, SymMatrixRef<T>, vector<double> > Solver_t;

  EigenSpaceDimension EigDim;
  //vector<T> Lambda;
  vector<double> Lambda;

  SymMatrixSolver():solver(NULL)  { }
  SymMatrixSolver(int n):solver(NULL) {EigDim.set(n);}
  SymMatrixSolver(int n, int m):solver(NULL) {EigDim.set(n,m);}
  ~SymMatrixSolver() {if(solver) delete solver;}

  //public interfaces for diagonalizations of a matrix
  void solve(T* ham, int ntot);
  void print(ostream& os) const;

private:
  
  Solver_t *solver;
};

template<class T>
void SymMatrixSolver<T>::solve(T* ham, int ntot) {

  if(ntot > 0) { // only when the dimention is non-zero

    EigDim.set(ntot,ntot); 

    if(!solver) solver = new Solver_t(EigDim);

    if(Lambda.size() < ntot) { // allocate eigen value vectors
      Lambda.erase(Lambda.begin(), Lambda.end());
      Lambda = vector<double>(ntot);
    }

    SymMatrixRef<T> matbase(ham);
    solver->getEigvectors(matbase, matbase, Lambda);
  }
}

template<class T>
void SymMatrixSolver<T>::print(ostream& os) const {
  os << "Eigen values of " << EigDim.N << " by " << EigDim.N 
     << " symmetric matrix:\n";
  for(int i=0; i<N; i++)
    os << setw(5) << i << setw(15) << setprecision(5) << Lambda[i] << endl;
}

#endif

/***************************************************************************
 * $RCSfile: SymMatrixSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: SymMatrixSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/
